Updated Spin and Orbital Parameters and Energy Dependent Pulse Behaviors of the Accreting Millisecond X-Ray Pulsar IGR J17591–2342

Soon after the discovery of millisecond pulsars from the radio band in the early 80s, low mass X-ray binaries (LMXBs) have been considered as their progenitors because a sufficient amount of angular momentum could transfer through the accretion disk to spin up the neutron stars and enable them to be millisecond pulsars during the evolution of LMXBs (Alpar et al. 1982; Radhakrishnan & Srinivasan 1982). Measuring the coherent pulsation in an LMXB provides an evolutionary link between radio millisecond pulsars and LMXBs. However, no coherent millisecond pulsation was detected in LMXBs until the discovery of the first accreting millisecond X-ray pulsar (AMXP), SAX J1808.4–3658, with a spin period of 2.5 ms (Chakrabarty & Morgan 1998; Wijnands & van der Klis 1998). To date, more than 20 LMXBs have been identified as AMXPs (Campana & Di Salvo 2018) with pulsation periods ranging from 1.67 ms (IGR J00291+5934, Galloway et al. 2005) to 9.51 ms (IGR J16597–3704, Sanna et al. 2018a).

The transient X-ray source IGR J17591–2342 was discovered by the INTernational Gamma-Ray Astrophysics Laboratory (INTEGRAL) on 2018 August 10 when it was scanning the Galactic center (Ducci et al. 2018). The follow-up observations were then made in X-ray (Bozzo et al. 2018; Ferrigno et al. 2018; Nowak et al. 2018), optical (Russell & Lewis 2018), near-infrared (Shaw et al. 2018), and radio bands (Russell et al. 2018). The Australia Telescope Compact Array (ATCA) made the observations on 2018 August 14, 19, and 25. In addition to reporting the precise source location, Russell et al. (2018) found that IGR J17591–2342 is a radio-loud system in comparison to other accreting neutron stars in their outburst states although no radio pulsation was detected. The Chandra made a 20 ks follow-up observation on 2018 August 23 (Nowak et al. 2018). The spectrum detected by the High Energy Transmission Gratings Spectrometer indicates that there is an outflow with a velocity of ∼2800 km s⁻¹ from the system and the discovery of the overabundance of calcium implies that IGR J17591–2342 could be formed from the collapse of a white dwarf system in a calcium-rich SN Ib (Nowak et al. 2019). Krimm et al. (2018) pointed out that the source had been detected by Swift/BAT on 2018 July 22, 20 days before its discovery by INTEGRAL and the flux reached its peak on July 25. Sanchez-Fernandez et al. (2018) and Ray et al. (2018) reported that the source experienced two rebrightenings starting on August 18 (MJD 58348) and September 6 (MJD 58366), respectively. Ferrigno et al. (2018) identified that IGR J17591–2342 is an AMXP with a pulsation frequency of 527 Hz (P_s ≃ 1.9 ms), an orbital period of 0.37 day (∼32,000 s), and a projected semimajor axis of 1.23 lt-s. More extensive studies were conducted by Sanna et al. (2018b) using the data collected by Nuclear Spectroscopic Telescope Array (NuSTAR) on 2018 August 12 and Neutron Star Interior Composition ExploreR (NICER) from 2018 August 15 to 24. In addition to deriving more precise spin and orbital parameters of the system, the broadband X-ray spectral analysis combining SWIFT and INTEGRAL data was also made.

The primary scientific objective of NICER is to investigate the interior composition of neutron stars by constraining their masses and radii using the precise measurement of the pulse profiles but AMXPs are not suitable for this goal due to their complex accretion flows (Gendreau et al. 2012). Nevertheless, the pulse profiles of the AMXPs can give us clues to understand the natures of the systems, such as the emission mechanism of the pulsations. Obtaining the precise pulse profile depends on accurate spin and orbital parameters for an AMXP. In this paper, we present our analysis results of timing and pulse profile properties of IGR J17591–2342 using the data collected by NICER of its 2018 outburst. The NICER observations and data reductions are briefly introduced in Section 2. In Section 3, we demonstrate the data analyses and results, including updating the spin and orbital parameters (Section 3.1), the energy dependent pulse behaviors (Section 3.2) and the spectral analysis (Section 3.3). Finally, the analysis results, including the implications of the neutron star magnetic field from the measured spin frequency derivative and the energy dependent pulse behaviors, are discussed in Section 4.

NICER started observing the IGR J17591–2342 on 2018 August 14 (ObsID:1200310101), soon after an outburst alert of the source was issued (Ducci et al. 2018), until near the end of the outburst on 2018 October 17 (ObsID:1200310139). The data files being analyzed in this work were first processed by the standard data filtering criteria described in the NICER Mission Guide³ with NICERDAS version 4.0 and the intervals with unexpected background flares were also excluded. Using the tool barycorr version 2.1 and by applying the JPL solar planetary ephemeris DE430, we corrected all the photon arrival times of events to the barycenter of the solar system according to the source position of the radio counterpart obtained by the ATCA observations (Russell et al. 2018). To reduce possible noise and background, only the events between energy range 1 to 12 keV (PI values between 100 and 1200) with both slow and fast chains triggered were selected for timing analysis. Figure 1 shows that the light curve of IGR J17591–2342 detected by NICER 1–12 keV band with the total exposure time of ∼91.5 ks. Two rebrightenings can be clearly seen in the light curve. More than the light curve shown in Ray et al. (2018), X-ray flux reaches its peak again around 2018 September 19–23 (MJD 58380–58384) in the second rebrightening and then decays close to zero in the following ∼20 days. No type-I burst was detected for all of the NICER observations.

Zoom In Zoom Out Reset image size

3.1. Timing Analysis

To update the spin and orbital parameters, all the selected events as described in Section 2 were divided into ∼300 s data segments for subsequent analysis. To choose the data segments with significant pulsation detection, we made the power spectra for all data segments using the ${Z}_{1}^{2}$ test (Buccheri et al. 1983) based on the circular orbital parameters proposed by Sanna et al. (2018b) derived from the NICER observations made between 2018 August 15 and 24. Only the data segments with a more than 99.9877% confidence level (${Z}_{1}^{2}\gt 18$) in pulsation detection were selected for further analysis. We found that all the data segments observed in August and September have significant pulsation detections but no pulsation can be found in the October observations, probably due to the very low count rates. In total, 264 data segments were selected for the following analysis.

To refine the spin and orbital parameters, a method similar to the one proposed by Sanna et al. (2016) was adopted. Suppose the orbital is circular, the relation between photon emission and arrival times, using the guess orbital parameters, can be written as

$$\begin{eqnarray}&&t^{\prime} =t-{A}^{(c)}\sin \left[2\pi {f}_{\mathrm{orb}}^{(c)}(t^{\prime} -{T}_{\mathrm{nod}}^{(c)})\right],\end{eqnarray} \tag{ 1 }$$

where $t^{\prime}$ is the photon emission time, t is the photon arrival time, $A\equiv {a}_{x}\sin i/c$ where ${a}_{x}\sin i$ is the projected orbital radius, f_orb is orbital frequency, T_nod is the time of passage of ascending node, and the parameters with superscript (c) represent the guess (calculated) parameters. The constant time delay between the barycenters of binary and solar system, D/c, is ignored in Equation (1). The photon emission time t' can be solved numerically using the iteration method described in Sanna et al. (2016).

For a constant spin frequency pulsar, the cycle count of the pulsar can be evaluated as

$$\begin{eqnarray}&&{N}_{c}={\nu }_{0}^{(c)}(t^{\prime} -{T}_{0}),\end{eqnarray} \tag{ 2 }$$

where ${\nu }_{0}^{(c)}$ is the guess spin frequency and T₀ is an arbitrary constant reference epoch. If there are small deviations between true and guess parameters, taking the first-order approximation, the pulse phase drift can be described as

$$\begin{eqnarray}\begin{array}{rcl}\delta \phi (t^{\prime} ) & \approx & -{\left.\left(\displaystyle \frac{\partial {N}_{c}}{\partial {\nu }_{0}}\right)\right|}_{(c)}{\left.\delta {\nu }_{0}-\left(\displaystyle \frac{\partial {N}_{c}}{\partial A}\right)\right|}_{(c)}\delta A\\ & & -{\left.\left(\displaystyle \frac{\partial {N}_{c}}{\partial {f}_{\mathrm{orb}}}\right)\right|}_{(c)}\delta {f}_{\mathrm{orb}}-{\left.\left(\displaystyle \frac{\partial {N}_{c}}{\partial {T}_{\mathrm{nod}}}\right)\right|}_{(c)}\delta {T}_{\mathrm{nod}},\end{array}\end{eqnarray} \tag{ 3 }$$

where δ are the differences between true and guess parameters, for example, $\delta {\nu }_{0}={\nu }_{0}-{\nu }_{0}^{(c)}$. For the orbital parameter parts in Equation (3), the evaluated photon emission time t' from Equation (1) is a function of guess orbital parameters, that is, $t^{\prime} =t^{\prime} ({A}^{(c)},{f}_{\mathrm{orb}}^{(c)},{T}_{\mathrm{nod}}^{(c)})$. Using Equation (2), the partial differential on the right-hand side of Equation (3), taking the parameter A as an example, should be

$$\begin{eqnarray}&&{\left.\left(\displaystyle \frac{\partial {N}_{c}}{\partial A}\right)\right|}_{(c)}={\nu }_{0}^{(c)}{\left.\left(\displaystyle \frac{\partial t^{\prime} }{\partial A}\right)\right|}_{(c)}\end{eqnarray} \tag{ 4 }$$

and from Equation (1)

$$\begin{eqnarray}&&\begin{array}{l}{\left.\left(\displaystyle \frac{\partial t^{\prime} }{\partial A}\right)\right|}_{(c)}=-\sin \left[2\pi {f}_{\mathrm{orb}}^{(c)}(t^{\prime} -{T}_{\mathrm{nod}}^{(c)})\right]\\ \quad -2\pi {f}_{\mathrm{orb}}^{(c)}{A}^{(c)}\cos \left[2\pi {f}_{\mathrm{orb}}^{(c)}(t^{\prime} -{T}_{\mathrm{nod}}^{(c)})\right]{\left.\left(\displaystyle \frac{\partial t^{\prime} }{\partial A}\right)\right|}_{(c)}.\end{array}\end{eqnarray} \tag{ 5 }$$

Therefore, solving ($\partial t^{\prime} /\partial A$) in Equation (5) and substituting it to Equation (4), we obtain

$$\begin{eqnarray}&&{\left.\left(\displaystyle \frac{\partial {N}_{c}}{\partial A}\right)\right|}_{(c)}=-\displaystyle \frac{{\nu }_{0}^{(c)}\sin \left[2\pi {f}_{\mathrm{orb}}^{(c)}(t^{\prime} -{T}_{\mathrm{nod}}^{(c)})\right]}{1+2\pi {f}_{\mathrm{orb}}^{(c)}{A}^{(c)}\cos \left[2\pi {f}_{\mathrm{orb}}^{(c)}(t^{\prime} -{T}_{\mathrm{nod}}^{(c)})\right]}.\end{eqnarray} \tag{ 6 }$$

Similar derivation can be applied to the other orbital parameters, f_orb and T_nod, and the pulse phase evolution can be expressed as

$$\begin{eqnarray}\begin{array}{rcl}\phi ({t}_{i}) & = & {\phi }_{0}-({t}_{i}-{T}_{0})\delta \nu \\ & & +\displaystyle \frac{{\nu }_{0}^{(c)}\sin \left[2\pi {f}_{\mathrm{orb}}^{(c)}({t}_{i}-{T}_{\mathrm{nod}}^{(c)})\right]}{1+2\pi {f}_{\mathrm{orb}}^{(c)}{A}^{(c)}\cos \left[2\pi {f}_{\mathrm{orb}}^{(c)}({t}_{i}-{T}_{\mathrm{nod}}^{(c)})\right]}\delta A\\ & & +\displaystyle \frac{2\pi {A}^{(c)}{\nu }_{0}^{(c)}({t}_{i}-{T}_{\mathrm{nod}}^{(c)})\cos \left[2\pi {f}_{\mathrm{orb}}^{(c)}({t}_{i}-{T}_{\mathrm{nod}}^{(c)})\right]}{1+2\pi {f}_{\mathrm{orb}}^{(c)}{A}^{(c)}\cos \left[2\pi {f}_{\mathrm{orb}}^{(c)}({t}_{i}-{T}_{\mathrm{nod}}^{(c)})\right]}\delta {f}_{\mathrm{orb}}\\ & & -\displaystyle \frac{2\pi {A}^{(c)}{f}_{\mathrm{orb}}^{(c)}{\nu }_{0}^{(c)}\cos \left[2\pi {f}_{\mathrm{orb}}^{(c)}({t}_{i}-{T}_{\mathrm{nod}}^{(c)})\right]}{1+2\pi {f}_{\mathrm{orb}}^{(c)}{A}^{(c)}\cos \left[2\pi {f}_{\mathrm{orb}}^{(c)}({t}_{i}-{T}_{\mathrm{nod}}^{(c)})\right]}\delta {T}_{\mathrm{nod}}\end{array}\end{eqnarray} \tag{ 7 }$$

where t_i is the time of ith data segment and ϕ₀ is a constant phase that depends on the reference epoch T₀. By fitting the evolution of pulse phases with Equation (7), the refined parameters can be derived by the corrections of the parameters. Taking the refined parameters as guess parameters, this correcting process can be iterated until the corrections of parameters become insignificant (e.g., smaller than the corresponding uncertainties of parameters).

Considering the spin frequency derivative, $\dot{\nu }$, Equation (2) should be rewritten as ${N}_{c}={\nu }_{0}^{(c)}(t^{\prime} -{T}_{0})+1/2{\dot{\nu }}^{(c)}{\left(t^{\prime} -{T}_{0}\right)}^{2}$ and the first-order approximation of the pulse phase drift should be written as

$$\begin{eqnarray}&&\begin{array}{l}\phi ({t}_{i})={\phi }_{0}-({t}_{i}-{T}_{0})\delta \nu -\displaystyle \frac{1}{2}{\left({t}_{i}-{T}_{0}\right)}^{2}\delta \dot{\nu }\\ +\,\displaystyle \frac{\left[{\nu }_{0}^{(c)}+{\dot{\nu }}^{(c)}({t}_{i}-{T}_{0})\right]\sin \left[2\pi {f}_{\mathrm{orb}}^{(c)}({t}_{i}-{T}_{\mathrm{nod}}^{(c)})\right]}{1+2\pi {f}_{\mathrm{orb}}^{(c)}{A}^{(c)}\cos \left[2\pi {f}_{\mathrm{orb}}^{(c)}({t}_{i}-{T}_{\mathrm{nod}}^{(c)})\right]}\delta A\\ +\,\displaystyle \frac{2\pi {A}^{(c)}\left[{\nu }_{0}^{(c)}+{\dot{\nu }}^{(c)}({t}_{i}-{T}_{0})\right](t^{\prime} -{T}_{\mathrm{nod}}^{(c)})\cos \left[2\pi {f}_{\mathrm{orb}}^{(c)}({t}_{i}-{T}_{\mathrm{nod}}^{(c)})\right]}{1+2\pi {f}_{\mathrm{orb}}^{(c)}{A}^{(c)}\cos \left[2\pi {f}_{\mathrm{orb}}^{(c)}({t}_{i}-{T}_{\mathrm{nod}}^{(c)})\right]}\delta {f}_{\mathrm{orb}}\\ -\,\displaystyle \frac{2\pi {A}^{(c)}{f}_{\mathrm{orb}}^{(c)}\left[{\nu }_{0}^{(c)}+{\dot{\nu }}^{(c)}({t}_{i}-{T}_{0})\right]\cos \left[2\pi {f}_{\mathrm{orb}}^{(c)}({t}_{i}-{T}_{\mathrm{nod}}^{(c)})\right]}{1+2\pi {f}_{\mathrm{orb}}^{(c)}{A}^{(c)}\cos \left[2\pi {f}_{\mathrm{orb}}^{(c)}({t}_{i}-{T}_{\mathrm{nod}}^{(c)})\right]}\delta {T}_{\mathrm{nod}}.\end{array}\end{eqnarray} \tag{ 8 }$$

To allow small eccentricity in the orbital model, two additional terms, $(e\cos \omega /2){A}^{(c)}[{\nu }_{0}^{(c)}+{\dot{\nu }}^{(c)}({t}_{i}-{T}_{0})]\sin [4\pi {f}_{\mathrm{orb}}^{(c)}({t}_{i}-{T}_{\mathrm{nod}}^{(c)})]$ and $-(e\sin \omega /2){A}^{(c)}[{\nu }_{0}^{(c)}+{\dot{\nu }}^{(c)}({t}_{i}-{T}_{0})]\cos [4\pi {f}_{\mathrm{orb}}^{(c)}({t}_{i}-{T}_{\mathrm{nod}}^{(c)})]$, are added into Equation (8) for the fitting, where e is eccentricity and ω is the longitude of the periastron.

The event times of each selected data segment were first folded with the initial spin and circular orbital parameters proposed by Sanna et al. (2018b), evaluated from 2018 August 15–24 NICER observations, and the phases were subsequently binned into 32 bins for a cycle to obtain the pulse profile. We found that all the pulse profiles can be well-fitted with a two-component sinusoidal function, that is $r(\phi )\,={a}_{0}+{\sum }_{k=1}^{2}[{a}_{k}\cos (2\pi k\phi )+{b}_{k}\sin (2\pi k\phi )]$. The peak of the best-fitted profile was selected as the fiducial point of the pulse phase for each data segment. The error of the pulse phase was evaluated using 10³ runs of the Monte Carlo simulation. If the initial parameters were correct, the pulse phase should have been consistent with being a constant within the measurement errors. The pulse phase evolution folded by the initial guess parameters is shown as Figure 2. A long-term trend of phase evolution can be clearly seen, that indicates that the spin and orbital parameters have to be further refined.

Zoom In Zoom Out Reset image size

The refining procedure as described above was adopted using the model containing spin (ν and $\dot{\nu }$) and orbital parameters (hereafter Model 1) to update the parameters by selecting T₀ on MJD 58338 that is close to the arrival time of the first event of NICER observations. Figure 2 shows the pulse phase residuals before and after completion of the refining process and the updated spin and orbital parameters obtained by NICER observations are listed in Table 1. The errors of spin parameters were quadratically added with the errors due to uncertainties of the source position proposed by Russell et al. (2018). A clear spin frequency derivative with a spin-down rate of $(-7\pm 1)\,\times {10}^{-14}$ Hz s⁻¹ was detected from NICER 2018 August to September observations. The eccentricity parts were further added into the model and found no significant eccentricity can be detected with a 2σ upper limit of 1 × 10⁻⁵.

Table 1.  Spin and Orbital Parameters of IGR J17591–2342 from NICER Observations

Parameter	Model 1a	Model 2b
Orbital period, Porb (s)	31684.7499(4)	31684.7499(4)
Projected semimajor axis, ${a}_{x}\sin i$ (lt-s)	1.227729(3)	1.227723(3)
Ascending node passage, Tnod (MJD/TDB)	58345.1719768(3)	58345.1719769(3)
Eccentricity, e	< 1 × 10−5	< 1 × 10−5
Spin frequency, ν0 (Hz)	527.42570061(2)	527.42570057(2)
Spin frequency derivative, $\dot{\nu }$ (Hz s−1)	$(-7\pm 1)\times {10}^{-14}$	$(-6\pm 1)\times {10}^{-14}$
Epoch of ν0, T0 (MJD/TDB)	58338.0	58338.0
α	⋯	0.73 ± 0.06
χ2/dof	626.94/258	475.45/257

Notes.

^aThe flux-dependent timing noise component is not included in the model. ^bThe flux-dependent timing noise component is included in the model.

Download table as:  ASCIITypeset image

However, the fitting from Model 1 gives a large reduced ${\chi }^{2}$ value of 2.43 for 258 degrees of freedom (dof) that is not statistically acceptable. Such a large deviation may be caused by the timing noise due to the movement of the hot spot on the neutron star surface that is highly correlated with X-ray flux observed in many of AMXPs (Patruno et al. 2009c). Furthermore, this phase drift may also induce a spurious spin frequency derivative in detection (Patruno et al. 2009c). To verify if this flux-dependent timing noise makes the large deviations, we first evaluated the linear Pearson correlation coefficient between the pulse phase residuals and the corresponding X-ray count rates. A correlation coefficient of −0.34 for 264 data points with a p-value of 1.3 × 10⁻⁸ was obtained, that implies the high correlation between the pulse phase residuals and the X-ray count rates (Figure 3). To minimize the effect from the flux-dependent timing noise, we adopted the numerical simulation results from Kulkarni & Romanova (2013) that the azimuthal location of the hot spot is proportional to ${\dot{M}}^{-1/5}$, where $\dot{M}$ is the accretion rate (also see Bult et al. 2020). Assuming that the 1–12 keV count rate is about proportional to the accretion rate, the phase drift due to the hot spot movement can be written as ${\rm{\Delta }}{\phi }_{n}=\alpha {r}^{-1/5}$ where r is count rate and α is a constant. We therefore added an additional term, $-\alpha {{r}_{i}}^{-1/5}$ where r_i is the count rate of the ith data segment, to Equation (8). The optimal spin and orbital parameters, as well as the proportional constant α, for this revised model (hereafter Model 2) are listed in Table 1. The spin frequency derivative $\dot{\nu }$ changes from $(-7\pm 1)\,\times {10}^{-14}$ Hz s⁻¹ to $(-6\pm 1)\times {10}^{-14}$ Hz s⁻¹, which implies the $\dot{\nu }$ evaluated from Model 1 is likely partly attributed to the flux-dependent timing noise. Considered together with this timing noise component, Figure 4 shows the quadratic trend of the pulse phase residuals folded by the constant spin frequency model plus the Keplerian orbit with the orbital parameters and α fixed at their optimal values listed in Table 1 (Model 2). Although the reduced χ² value decreases significantly from 2.43 to 1.85, that indicates the revised model is better than the previous one, the value is still too large to be statistically acceptable. However, there is no structure in the residuals. This poor fitting is likely due to the remaining unmodeled timing noise (Patruno & Watts 2012).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

On the other hand, we cannot exclude the possibility that the quadratic trend in the pulse phase evolution is due to red noise. To verify this, a simulation of a first-order autoregressive (AR(1)) process as ${\phi }_{i}={a}_{i}{\phi }_{i-1}+{\varepsilon }_{i}$ was made. First, assuming that the quadratic trend on Figure 4 (with $\dot{\nu }=-6\times {10}^{-14}$ Hz s⁻¹) is purely due to the AR(1) process, we selected the equal-spacing points from the curve and evaluated a_i = a = constant using the autocorrelation function. Furthermore, suppose that the process is damped random walk, $a=\exp (-{\rm{\Delta }}t/\tau )$, where Δt is the time interval of the two consecutive equal-spacing points. The damping timescale τ can be obtained. Applying the τ to the unevenly spacing data points of pulse phase evolution, we simulated an AR(1) process as ${\phi }_{i}\,=\exp [-({t}_{i}-{t}_{i-1})/\tau ]{\phi }_{i-1}+{\varepsilon }_{i}$, where t_i is the time of the ith data segment and ${\varepsilon }_{i}$ is pure Gaussian white noise with mean zero and the variance that equals the square of phase error of the ith data segment, and then fitted a quadratic curve to the simulated data to find the probability that $| \dot{\nu }| \gt 6\,\times {10}^{-14}$ Hz s⁻¹ and ${\chi }_{\nu }^{2}\lt 1.85$. After 10⁶ runs of the Monte Carlo simulation, none of them is satisfied by the criteria above. We therefore conclude that the quadratic trend is unlikely caused by the red noise with a probability less than 10⁻⁶.

3.2. Energy Dependent Pulse Behaviors

The combined pulse profile of all the data segments with significant pulse detections folded with the optimal spin and orbit parameters listed in Table 1 is shown as Figure 5 with an rms pulsed fractional amplitude of 8.8%. Fitting the multiple sinusoidal functions to the pulse profile shows that it needs a four-component sinusoidal function, that is, $r(\phi )={a}_{0}+{\sum }_{k=1}^{4}[{a}_{k}\cos (2\pi k\phi )+{b}_{k}\sin (2\pi k\phi )]$, to obtain an acceptable fitting with ${\chi }_{\nu }^{2}=0.98$. The rms pulsed fractional amplitude of the four components are 7.95%, 3.66%, 0.60%, and 0.18% for the fundamental, second, third, and fourth harmonics, respectively.

Zoom In Zoom Out Reset image size

To investigate the energy dependent pulse behaviors, the 1–12 keV events were divided into 11 energy bands. Each of them has roughly the same event numbers, except for the bands 10 and 11 that have only about half of the event numbers compared with the other bands for further investigation of the pulse behaviors of the harder X-ray bands. The pulse profiles of these 11 bands are shown in Figure 6. All pulse profiles are well-fitted with a three-component sinusoidal function. To study the energy dependent pulse arrival time delay, that has been seen in many other AMXPs (Patruno & Watts 2012), we applied the cross-correlation between the best-fitted pulse profiles and that of the softest band (band 1, 1–1.64 keV). The uncertainties of the delays were evaluated by 10³ runs of the Monte Carlo simulation, except for band 1, which directly used the uncertainty of its fiducial point. Figure 7 shows the energy dependent pulse arrival delay relative to band 1. There is no significant soft lag up to band 7 (3.74 keV) and pulsations of harder energy bands lead the band 1 up to ∼0.06 cycles (∼114 μs) at the hardest band (band 11, 7.51–12 keV). There is neither soft lag saturation (e.g., HETE J1900.1–2455, Galloway et al. 2007) nor hard lags (e.g., IGR J00291+5934, Falanga et al. 2005) that can be seen in IGR J17591–2342 from the NICER observations. Figure 8 shows the rms pulsed fractional amplitudes of different energy bands. They first monotonically increase from ∼4.7% for band 1 up to ∼11.5% for band 9 (4.41–5.55 keV) and then decrease to ∼5.9% in the hardest energy band (band 11, 7.51–12 keV). This implies that there may be soft unpulsed emissions in the spectrum that will be further discussed in Sections 3.3 and 4.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.3. Energy Spectrum

The energy dependent pulse behaviors of an AMXP would be highly related to its energy spectrum (see Section 4 for more detailed discussions), so we performed spectral analysis of the data collected by NICER. The spectrum of energy range from 0.7 to 12 keV was extracted from the same data set as the one used for timing analysis. The background was evaluated from NICER observations of Rossi X-ray Timing Explorer blank field region 5 (Jahoda et al. 2006), followed by the same filtering as we processed to the source data. To explain the energy dependent pulse arrival time delay detected in SAX J1808.4–3658, Gierliński et al. (2002) and Poutanen & Gierliński (2003) proposed a two-component model composed of blackbody emission from a hot spot and Comptonized emission from a hot slab, to describe the pulsed emission from an AXMP. Using XSPEC version 12.10 (Arnaud 1996), the spectrum was first modeled with an absorbed blackbody plus a Comptonized continuum as tbabs(bbodyrad+nthcomp). For the Comptonized component, we found it is hard to well constrain the electron temperature, which is likely due to the lower energy range (0.7–12 keV) for the observed spectrum. The electron temperature was therefore fixed at 22 keV, same as the one found in IGR J17591–2342 by Sanna et al. (2018b) from the broadband spectral fitting. However, a large reduced χ² value of 3.13 (dof = 1125) was yielded from the fitting. A majority deviation between the data and the model was found in the low energy range (<1.4 keV).

The energy dependent pulsed fractional amplitudes (Figure 8) imply that there could be a soft, unpulsed component in the spectrum. Therefore, a multicolor disk blackbody component was added into the spectral model, that is, tbabs(diskbb+bbodyrad+nthcomp). Adding this additional component significantly improved the fitting with a reduced χ² of 1.54 (dof = 1123) and the spectral parameters are presented in Table 2. Figure 9 shows the unfolded spectrum of the NICER data set we analyzed. Although this reduced χ² is still statistically unacceptable, the large change of the χ² (Δ χ² = 1792) implies that the unplused disk blackbody is an essential component for the spectral model. On the other hand, we noticed that the remaining residuals are mostly concentrated below ∼2.5 keV, that is likely caused by the instrumental residuals in the NICER response function. In addition to the disk blackbody component, it is evident that the emission from the blackbody component dominates the one from the Comptonized component for photon energies ≲3.2 keV. This characteristic is related to the energy dependent pulse arrival time that will be further discussed in Section 4.

Zoom In Zoom Out Reset image size

Table 2.  The Best-fitted Spectral Parameters of IGR J17591–2342 from NICER Observations

Parameter	tbabs(diskbb+bbodyrad+nthcomp)
NH (1022 cm−2)	${2.01}_{-0.02}^{+0.03}$
kTin (keV)	0.081 ± 0.002
Disk flux (1–12 keV) (erg s−1 cm−2)	${9.20}_{-0.87}^{+1.04}\times {10}^{-12}$
kTBB (keV)	0.64 ± 0.02
NormBBa	${57.7}_{-3.0}^{+3.5}$
BB flux (1–12 keV) (erg s−1 cm−2)	${9.42}_{-0.83}^{+0.64}\times {10}^{-11}$
BB flux (1–4 keV) (erg s−1 cm−2)	${8.14}_{-0.61}^{+0.45}\times {10}^{-11}$
Γ	${2.85}_{-0.50}^{+0.92}$
kTe(keV)	22 (fixed)
kTseed(keV)	1.48 ± 0.19
nthcomp flux (1–12 keV) (erg s−1 cm−2)	${2.15}_{-0.14}^{+0.08}\times {10}^{-10}$
nthcomp flux (1–4 keV) (erg s−1 cm−2)	${5.13}_{-0.49}^{+0.76}\times {10}^{-11}$
${\chi }_{\nu }^{2}$ (dof)	1.54 (1123)

Note.

^aThe normalization of bbodyrad is presented in units of ${R}_{\mathrm{km}}^{2}/{D}_{10\,\mathrm{kpc}}^{2}$, where R²_km is blackbody emission radius in the units of km and ${D}_{10\,\mathrm{kpc}}^{2}$ is the source distance in the unit of 10 kpc.

Download table as:  ASCIITypeset image

We have updated the spin and orbital parameters of AMXP IGR J17591–2342 from the data collected by NICER in 2018 August and September. The updated parameters are generally consistent with those reported by Sanna et al. (2018b) but with more precise values (see Table 1). The only exception is the spin frequency derivative. We found that the neutron star in IGR J17591–2342 is, averagely speaking, spun-down with a rate of −6 × 10⁻¹⁴ Hz s⁻¹ during the entire duration of NICER observations, instead of being a constant spin frequency (with an insignificant spin-up rate of $(2.1\pm 1.6)\times {10}^{-13}$ Hz s⁻¹) as proposed by Sanna et al. (2018b) from 2018 August 15 to 24 observations. This value essentially agrees with the theoretical prediction (∼10⁻¹³ Hz s⁻¹). The spin frequency derivative allows us to roughly estimate the magnetic field on the surface of a neutron star in this AMXP. The neutron star can be spun-up by accretion and be spun-down due to the magnetic drag in the accretion disk and the magnetic dipole radiation. Using Equation (23) in Rappaport et al. (2004) for fast X-ray pulsars but also considering the magnetic dipole radiation, the torque acting on the neutron star can be rewritten as

$$\begin{eqnarray}&&{\rm{\Gamma }}=2\pi I\dot{\nu }\simeq \dot{M}\sqrt{{{GMr}}_{c}}-\displaystyle \frac{{\mu }^{2}}{9{r}_{c}^{3}}-\displaystyle \frac{2}{3}{\mu }^{2}{\sin }^{2}\alpha {\left(\displaystyle \frac{2\pi \nu }{c}\right)}^{3},\end{eqnarray} \tag{ 9 }$$

where I is the moment of inertia of the neutron star, $\dot{M}$ is the mass accretion rate, M is the neutron start mass, μ is the magnetic dipole moment, r_c is the corotation radius that the radius with the Keplerian frequency is equal to the spin frequency, α is the angle between rotational axis and magnetic dipole, G is the gravitational constant and c is the speed of light. Using the numerical values of $\dot{M}\simeq 5.2\times {10}^{-10}{M}_{\odot }\,{\mathrm{yr}}^{-1}$ evaluated from the broadband spectrum (Sanna et al. 2018b), M = 1.4M_⊙, ${r}_{c}={\left({GM}/4{\pi }^{2}{\nu }^{2}\right)}^{1/3}$ = 25.7 km, I ≃ 10⁴⁵ g cm², and ν and $\dot{\nu }$ from Table 1, the magnetic dipole moment could be estimated as $\mu \simeq (3.8\mbox{--}4.1)\times {10}^{26}$ G cm³ for $\sin \alpha$ = 1 to 0, respectively, which implies that the magnetic field B ≃ 4 × 10⁸ G on the neutron star for a radius of 10 km, consistent with the magnetic field of 10⁸–10⁹ G for millisecond pulsars. However, we note that the spin frequency derivative contributed by the magnetic dipole radiation can be as high as −2.2 × 10⁻¹⁴ Hz s⁻¹ (for $\sin \alpha =1$), ∼37% of the total spin frequency derivative. If we assume that the magnetospheric radius ${r}_{m}\simeq {r}_{c}=\xi {r}_{{\rm{A}}}$, where ${r}_{{\rm{A}}}={\left({\mu }^{4}/2{GM}{\dot{M}}^{2}]\right)}^{1/7}$ is the Alfvén radius, forms the numerical values above, r_A ≃ 48.6 km was obtained and it indicates ξ ≃ 0.53, close to the typical value for disk accretion (ξ ∼ 0.5).

The energy dependent pulse arrival time can usually be observed in AMXPs. Soft lags are most often seen in AMXPs, where the pulse arrival times of the pulsations from softer energy bands lag compared to the ones from the harder energy bands. This phenomenon was first discovered by Cui et al. (1998) during the 1998 outburst of the SAX J1808.4–3658. The soft lag tendency can clearly be detected from 2 keV extended up to 10 keV (hereafter, the break point) for about 200 μs (∼0.08 cycles) but saturated for harder X-ray bands. Similar energy dependent pulse arrival time behavior can also be found in XTE J1751–305 (Gierliński & Poutanen 2005), XTE J1814–338 (Watts & Strohmayer 2006), HETE J1900.1–2455 (Galloway et al. 2007), and XTE J1807–294 (Chou et al. 2008) with different break point energies. However, the energy dependent pulse arrival time acts differently in IGR J00291+5934. Falanga et al. (2005) found that the soft lags can extend to the break point energy at ∼6 keV but the tendency reverses instead of saturation for the pulse of energy bands harder than the break point (hereafter hard lags) during its 2005 outburst. Similar behaviors are also observed during its 2015 outburst (Sanna et al. 2017) with a different break point energy at ∼8 keV. Such hard lags are also marginally detected in IGR J17511–3057 (Falanga et al. 2011) and IGR J17498–2921 (Falanga et al. 2012). On the other hand, the energy dependent pulse arrival times for two AMXPs are rather peculiar. The pulse arrival times seem independent of energy bands below ∼17 keV for SAX J1748.9–2021 (Patruno et al. 2009a) and only hard lags are detected in IGR J18245–2452 (De Falco et al. 2017). The newly discovered AMXP IGR J17379-3747 also shows unusual energy dependent pulse arrival times. Regular soft lags can be observed between 0.5 and 6 keV with no break point being detected; however, the phase difference between the softest and hardest bands can be as high as ∼0.3 cycles (or ∼110°; Bult et al. 2019).

A two-component model was proposed by Gierliński et al. (2002) and Poutanen & Gierliński (2003) to explain the soft lags seen in SAX J1808.4–3658 during its 1998 outburst. There are two emission components, a soft blackbody component from a hot spot on the surface of neutron star around the accretion channel and a hard Compotonized component from up-scattering of the cooler seed photons by the hot electrons in a slab of accretion shock above the hot spot. Due to the scattering, the angular distribution of Comptonized photons is wider than that of the blackbody. Considered with Doppler boostering from the fast rotating AMXP, and the light bending in the Schwarzchild geometry, the pulsation from the Comptonized component leads to the pulsation from the blackbody component. For the softest energy band, the blackbody component dominates the emissions. As the energies of photons increased, the Comptonized emissions become stronger and the blackbody component becomes weaker. This causes the pulsation in harder energy bands to precede those in soft (i.e., soft lags). Beyond the break point, the effect from the blackbody component can be completely neglected, leaving only the Comptonized component, so the energy dependent pulse arrival time reaches saturation. This two-component model can explain the energy dependent pulse arrival times of SAX J1808.4–3658 (Cui et al. 1998), XTE J1751–305 (Gierliński & Poutanen 2005), XTE J1814–338 (Watts & Strohmayer 2006), HETE J1900.1–2455 (Galloway et al. 2007), and XTE J1807–294 (Chou et al. 2008) but cannot be applied to hard lags seen in IGR J00291+5934 (Falanga et al. 2005; Sanna et al. 2017).

We attempt to adopt the two-component model to explain the energy dependent behaviors of pulsations presented in Section 3. Figure 7 shows that the pulse phases are nearly constant between 1 and ∼4 keV and the soft lags appear for the harder energy bands. This indicates that the blackbody emissions from the hot spot dominate the Comptonized component's emission for the photon energy below ∼4 keV. The spectral analysis results also support this implication. The blackbody component is stronger than the Comptonized component for photon energy ≲4 keV as shown in Figure 9. The blackbody flux is a factor of 1.6 (=8.14/5.13) larger than the Comptonized flux for the 1–4 keV band even though the Comptonized flux is a factor of 2.3 (=21.5/9.42) larger than the blackbody flux for the 1–12 keV band (see Table 2). This indicates that the IGR J17591–2342 has relatively strong blackbody emissions. A similar phenomenon can be observed in HETE J1900.1–2455 (Galloway et al. 2007, Figure 3), IGR J17511–3057 (Falanga et al. 2011, Figure 6), and IGR J17498–2921 (Falanga et al. 2012, Figure 5). By contrast, some AMXPs show sharp soft lag tendency starting at ∼2 keV, such as SAX J1808.4–3658 in the 1998 outburst (Cui et al. 1998, Figure 4) and IGR J00291+5934 in 2005 outburst (Falanga et al. 2005, Figure 6), which implies that these sources may have relatively weaker blackbody emissions from the hot spot.

We did not detect the break point in IGR J17591–2342 from its 2018 outburst between the 1 and 12 keV energy bands. This implies that the break point is likely ≳10 keV. From the two-component model, the break point energy indicates that the contribution from the blackbody component can be ignored beyond this energy. Some AMXPs have softer break point energies, such as XTE J1814–338 (at ∼6 keV, see Figure 2 in Watts & Strohmayer 2006), HETE J1900.1–2455 (at ∼6 keV, see Figure 3 in Galloway et al. 2007), IGR J00291+5934 in the 2005 outburst (at ∼6 keV, see Figure 6 in Falanga et al. 2005) and 2015 outburst (at ∼8 keV, see Figure 4 in Sanna et al. 2017), whereas a few AMXPs have harder break point energies, such as XTE J1751–305 (>10 keV, see Figure 8 in Gierliński & Poutanen 2005) and IGR J17511–3057 (>10 keV, see Figure 6 in Falanga et al. 2011). Falanga et al. (2011) argued that the harder break point energy could be due to higher blackbody component temperature. However, a similar effect could also be the result of a relatively strong blackbody emission that extends the contribution from the pulsation of the blackbody component to a harder energy band, such as what we observed in IGR J17591–2342.

On the other hand, the energy dependent pulse amplitude analysis shows that the rms pulsed fractional amplitudes increase from 1 (4.7%) to ∼5 keV (11.5%) and then decrease for higher energy bands (Figure 8). The decrease of pulse fractional amplitudes for softer energy bands below ∼5 keV is very likely due to the unpulsed component (i.e., the DC term) in the pulse profiles being stronger for the softer energy bands. Adding a multicolor disk blackbody component can significantly improve the spectral fitting, especially for the low energy part. This implies that the unpulsed component could mainly be due to the thermal emission from the accretion disk. This phenomenon is analogous to the one found in SAX J1808.4–3658 in its 2008 outburst from XMM-Newton observations (Patruno et al. 2009b). The pulsed fractional amplitudes largely decrease below 2 keV and a disk blackbody component is also required for acceptable spectral fitting. In fact, a similar phenomenon can also be observed in IGR J00291+5934 during its 2015 outburst for energy bands less than ∼2 keV (see Figure 4 in Sanna et al. 2017), although it seems that a disk blackbody component is unnecessary for their spectral fitting. However, we note that the effect of unpulsed emission can be extend up to ∼5 keV from the energy dependent pulsed fractional amplitudes but the disk blackbody component from the spectral fitting seems too small to have such a large influence on the unpulsed component of the pulse profile. This is probably due to the imperfect modeling of the Compotonized component, because the electron temperature is fixed at 22 keV. The other possibility is that there is other unpulsed emission from the system, such as thermal emission from the surface of neutron star in the region other than the hot spot. We have attempted to add an additional blackbody component to model the thermal emission from the other part of the surface of neutron star with smaller temperature than the hot spot but failed because these two blackbody components tend to merge together during the fittings, probably due to insufficient spectral resolution to distinguish these two components.

The authors thank C.-P. Hu for helpful discussions about the red noise. The NICER data for this research were obtained from the High Energy Astrophysics Science Archive Research Center (HEASARC) online service, provided by the NASA/Goddard Space Flight Center. This work was supported by the Ministry of Science and Technology of Taiwan through the grants MOST 108-2112-M-008-005- and MOST 109-2112-M-008-004-.

Facilities: ADS - , HEASARC - , NICER. -

Software: heasoft (HEASARC 2014), nicerdas (v004), XSPEC (Arnaud 1996).